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We introduce a new template for the detection of gravitational waves from compact binary systems 
which is based on Chebyshev polynomials of the first kind. As well as having excellent convergence 
properties, these polynomials are also very closely related to the elusive minimax polynomial. In 
this study we have limited ourselves to the test-mass regime, where we model a test particle in a 
circular equatorial orbit around a Schwarzschild black hole. Our objective is to model the numeri- 
cal gravitational wave flux function, starting with the post-Newtonian expansion from Black Hole 
If^ • Perturbation Theory. We introduce a new Chebyshev approximation to the flux function, which 

' due to a process called Chebyshev economization gives a better model than either post-Newtonian 

, or Pade based methods. A graphical examination of the new flux function shows that it gives an 

' excellent fit to the numerical fiux, but more importantly we find that at the last stable orbit the 

error is reduced, < 1.8%, at all orders of approximation. We also find that the templates constructed 
O , using the Chebyshev approximation give better fitting factors, in general > 0.99, and smaller errors, 

■ < 1/10%, in the estimation of the Chirp mass when compared to a fiducial exact waveform, con- 

structed using the numerical fiux and the exact expression for the orbital energy function, again at 
all orders of approximation. We also show that in the test-mass case, the new Chebyshev template is 
superior to both PN and Pade approximant templates, especially at lower orders of approximation. 
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I. INTRODUCTION 



(N 

The inspiral of compact objects such as neutron stars and black holes are expected to be a major source of 
^ ' gravitational waves (GW) for the ground-based detectors LIGO, VIRGO, GEO600 and TAMA [lllllli, as well as 

■ for the future planned space-based detector LISA Q . Due to the effect of radiation reaction, the orbit of the binary 
f-~>) I system slowly decays over time. As this happens the amplitude and frequency of the waveform increases emitting a 

«. 'chirp' waveform. 

The most prominently used method for detecting these waveforms is matched filtering 0. In this method one 
I . creates a series of theoretical waveforms or templates which depend on the parameters of the system. These templates 
are then cross correlated with the detector output, weighted with the inverse of the power spectral density. We 
. . , demand that the templates achieve a predetermined threshold in the root-mean-square correlation before having any 
confidence in the detection of a GW and the parameters it estimates. One of the main problems with matched filtering 
is that it is highly sensitive to the phase of the waveform. A template that goes out of phase with a possible signal, 
even for a small number of cycles, leads to a severe degradation in the recovered signal-to-noise ratio (SNR). It is 

■ therefore imperative that we know the phase of the wave to very good accuracy. 
There have been many efforts to create templates which will approximate a possible signal to high accuracy. 

On one hand we have the post-Newtonian (PN) expansion of Einstein's equations to treat the dynamics of the 
system !H S ^3 0; ^1 ■ This works well in the adiabatic or slow- motion approximation for all mass ranges. On 
the other hand we have black hole perturbation theory JA, 15^ il^ dlt tl^ JM which works for any velocity, but only 
in situations where the mass of one body is much greater than the other. While templates have been generated to 5.5 
PN order for a test-mass orbiting a Schwarzschild black hole and to 3.5 PN order for non-spinning binaries of 
comparable mass 12, 13J, a number of difhculties still need to be tackled. The main problem is that both templates 
are a function of the orbital energy and GW flux functions. In the test-mass case, an exact expression in known for 
the orbital energy, but we have a PN expansion for the flux function. In the comparable-mass case, a PN expansion is 
known for both functions. It has been shown that the convergence of both methods is too slow to be useful in creating 
templates that can be confidently used in a GW search [l^ |23, lU, 122. . We also know that the PN approximation 
begins to break down when the orbital separation of the two bodies is r < lOM This means that as we approach 
the Last Stable Orbit (LSO) the templates begin to go out of phase with a possible signal due to the increase of 
relativistic effects. It is because of this that we need to have the most accurate templates that we can possibly have, 
without being computationally expensive. 

As well as the PN templates, there are many other templates on the market at present, all of whom have certain 



advantages (and disadvantages) over others. It was shown [23, |25|, |26i| that templates based on resummation methods 
such as Pade approximation have a faster convergence in modelling the gravitational waveform. This gave a superior 
template in the test-mass and comparable-mass case (if the comparable-mass case can be treated as a smooth de- 
formation of the test-mass case. See, for example, Blanchet for arguments against). The Pade based templates 
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were then used to partially construct Effective One Body templates |29| which went beyond the adiabatic ap- 
proximation and modelled the waveform into the merger phase. Other more phenomenological templates such as the 
BCV 1^ 1^ templates seem to be excellent at detecting GW, but are not necessarily the best template to use 
in the extraction of parameters. In any ideal detection strategy we would like to have a single template which could 
provide a confident detection, while extracting the parameters of the system to very high accuracy. 



A. Improving Template Construction Using the Minimax and Chebyshev Polynomials. 

The question is can we improve on the construction of the templates. The answer, theoretically, is yes. The problem 
with expansions like a Taylor series is that they are based on Weierstrass's theorem. This states : For any continuous 
function f{x) G C[a,6] and for any given e > 0, there exists a polynomial p„{x) for some sufficiently large n such 
that \ f{x) — Pnix)\max < £• So as long as we expand a series to a sufficiently large number of terms, we can always 
approximate a function very well. However, this is not always possible as the number of necessary terms may be too 
high (for example, it takes about 5000 terms in an expansion of arctan{x) to deliver five significant figures at x equal 
to unity js^*). or we may be dealing with an approximation to some function where it is very difficult to calculate more 
than a few terms (e.g. the PN expansion of the flux function). We know that the flux function has been approximated 
by an expansion to 11 terms for a test-mass body, and 7 for a comparable mass body. We also know from previous 
studies that this may not be sufficient. A more promising possibility is based on the Chebyshev Alternation theorem 
for polynomials which states ''For any continuous function f{x) € C[a,b], a unique minimax polynomial approximation 
Pn{x) exists, and is uniquely characterised by the 'alternation property' (or 'equios dilation property') that there are 
(at least) n + 2 points in [a,b] at which \ f{x) — Pn{x)\ attains its maximum absolute value with alternating signs'. 
Thus, the reason the minimax polynomial is so sought after is that the Chebyshev Alternation theorem guarantees 
that there is only one such polynomial and it has the necessary condition of having an 'equal-ripple' error curve. 
Unfortunately, the minimax polynomial is extremely difficult, if not impossible, to find. 

However, all is not lost. The family of Ultraspherical (or Gegenbauer) polynomials are defined by 

P^,-\x) = C„ (1 - x^)-" ^ (1 - (-1 < a < 00) , (1) 

where C„ is a constant. These polynomials are orthogonal over x G [—1, 1] with respect to the weight function 

(1 — a;^)". A feature of the polynomials Pii^\x) is that they have n distinct and real zeros and exhibit an oscillatory 
behaviour in the interval [—1, 1]. For values of —1/2 < a < 00, the amplitude of the oscillations increase as we move 
from the origin to the endpoints. For —1 < a < —1/2, the oscillations decrease as we move to the endpoints. However, 
for a — —1/2 the amplitude of the oscillations remain constant throughout the interval. 

This value of a corresponds to the Chebyshev polynomials of the first kind, Tn{x), (hereafter Chebyshev polyno- 
mials) with C„ — (— l)"2"ri!/(2n)!. These polynomials are closely related to the minimax polynomial due to the 
fact that there are n -I- 1 points in [-1,1] where Tn{x) attains a maximum absolute value with alternating signs, i.e. 
|T„(a:)| = ±1 |34j . This is their main strength over the Taylor or Fade type expansions as neither obey the alternation 
theorem and are thus nowhere close to the minimax polynomial. In fact it can be shown (35| that the Chebyshev 
polynomials exhibit the fastest convergence properties of all of the Ultraspherical polynomials. The n -t- 1 extrema 



are given by 

a;fc=cos(— j fc = 0, ..,n, (2) 

with n zeros given by 

/[/c- i]7r\ 

Xk = cos fc = 1, .., n. (3) 

\ n J 

The Chebyshev polynomials are formally defined by 

Tn{x) = cos(n0) when x = cos(6'). (4) 

From de Moivre's theorem, we know that cos {nO) is a polynomial of degree n in cos {9), e.g. 

cos(Oe') = l, cos (6*) = cos (6*) , cos (26*) = 2cos^ (6*) - 1, , (5) 

which allows us to write the Chebyshev polynomials 

To{x) = 1, ri(x) = X, T2{x) - 2x^ - 1, (6) 
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Therefore, a Chebyshev series in x corresponds to a Fourier series in 6. This close relation to a trigonometric function 
means these polynomials are extremely useful in approximating functions. The Chebyshev polynomials are orthogonal 
polynomials with respect to the weight (l — x^) according to 
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T,{x)Tj{x) 
1 Vl~a;2 I ^ i=j = Q 



mx)\T,{x))= rf:. 7 f ^-jVO (7) 



with initial conditions 

Toix) = l , Ti{x)=x. (8) 

We calculate the higher order Chebyshev polynomials using the recurrance relation 

T„(x) = 2xTn-i{x) - Tn-2{x), (9) 

to calculate the higher order polynomials, given the initial conditions. 

For our purposes, we need to approximate polynomials which are a function of the dimensionless velocity v in the 
domain v g [0,Uiso], where viso = I/'s/B is the velocity at the LSO for a test-particle orbiting a Schwarzschild black 
hole. In this case we use the Shifted Chebyshev polynomials, designated T*{v). We can transform from the interval 
[1, —1] to an arbitrary interval [a, b] using 

^^2x_^a^ V:Ee [a,6],se hl,l]. (10) 

In this case we have 

s = — -l = V24i;-l V-y e [O,wj,o]. (11) 

We can now write the shifted Chebyshev polynomials in the form 

r„*(«) = r„(s) = Tn{V2iv - 1), (12) 

and the recurrence relation as 

T:{v) = 2(V2iv - \)T*_^{v) - T:_,{v), (13) 
such that the shifted polynomials are have the initial conditions 

T*iv)^l, T*{v)^V24v-l. (14) 
In Fi gure (yT|) we have plotted the first few Chebyshev polynomials in the interval v G [0, ^/so]- 



B. Organization of the paper. 

In Section^ we define the restricted FN waveform. While not being completely correct, this waveform suffices to 
allow us to test the modelling of the gravitational wave flux function. In this section we define the flux function and 
the orbital energy function for a test particle in an equatorial circular orbit around a Schwarzschild black hole. In 
Section IIIII we outline the steps involved in defining the Chebyshev approximation to the gravitational flux function 
to 5.5 FN order. We also show in this section how the monomials in the velocity v can be written in terms of 
Chebyshev polynomials. In Section llVI we first of all define Chebyshev Economization and examine the truncation 
errors involved in both the FN and Chebyshev approximations to the flux. We then look at the graphical modelling 
of the flux function and analyse how the new approximation competes with both the FN and Fade approximations 
when compared to a numerical flux. In particular we look at how the new approximation improves the fit to the 
numerical flux at the LSO. In Section we look at fitting factors and parameter estimation for a choice of three 
test-mass systems. In this section we again compare the performance of the Chebyshev, FN and Fade templates when 
compared to an 'exact' signal constructed using the exact orbital energy and a numerically calculated flux. 

Throughout the paper we use G = c = 1. 
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II. THE TEST-MASS GRAVITATIONAL WAVEFORM. 



In the TT gauge, the response of a detector to an impending gravitational wave is given by h{t) — h+F^ + h-^F^ 
where /i+,x are the polarizations of the wave and F^''^ are the beam pattern functions of the detector. For ground 
based detectors, the duration of the wave is so short (less than a minute in most cases), we can take the beam pattern 
functions to be constant. In this case we can write the restricted PN waveform as 



h{i) 



Arm] 



Cw2(i)cos($(t) + $o), 



(15) 



where m = mi + TO2 is the total mass of the system, 77 ~ 71111112 I'm? is the reduced mass ratio, d is the distance to the 
source, C is a constant, v{t) is the orbital velocity of the bodies, <&(<) is the phase of the wave and $0 can be taken as 
being the phase at the time of arrival of the wave in the detector. In the stationary phase approximation the Fourier 
transform for positive frequencies reads [s^ Isvi Is^ . fsflj l 



HI) 



h{t) exp(27ri/i) dt = AJ 



and, since h(t) is real, h{—f) = h*{f). Here the amplitude A is found by normalizing the waveform such that, 

df 



{h\h) = 2 



Snif) 



Hf)h*{f) + h*{f)h{f) 



= 1, 



(16) 



(17) 



where the * denotes complex conjugate and h{f) is the Fourier transform of h{t). Also, Sh{f) is the one-sided 
noise power spectral density of the detector. In this study we will use one of the projected third-generation detectors, 
EURO [i^l . The reason we use this particular curve is due to the fact that this detector has a predicted lower frequency 
cutoff of about 10 Hz allowing us to investigate larger mass systems that would have already begun coalescence in the 
bandwith of the first generation detectors. We have plotted the effective noise h = \/ fShif) for the main first and 
third generation detectors in Figure Q. Solving the above equation for A gives 



A^ 



df 



-7/3 



-1/2 



The phase of the Fourier transform in the stationary phase approximation, ■(/'(/) i is given by 



E'jv) 
F{v) 



dv. 



(18) 



(19) 
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where (to, (j)a) can be taken as being the time of arrival and the phase of the wave at time of arrival of the gravitational 
wave, Vf — (Trm/)^/'^ is the instantaneous velocity, E'{v) = dE/dv is the derivative of the orbital energy with respect 
to the velocity v = (mfi)^/"^ = x^^"^, where fl is the angular velocity as observed at infinity and x is an invariant 
velocity parameter observed at infinity. Finally, F(v) is the gravitational wave flux function. Instead of solving the 
integrals in the above equation it is numerically more efficient to solve the following equivalent differential equations 
for the phasing formula in the Fourier domain: 



2TTt = 0. 



dt 



ttA'P E'if) 



0. 



(20) 



For a test-mass particle in circular orbit about a Schwarzschild black hole, an exact expression for the orbital energy 
exists For our purposes we need it's derivative with respect to the velocity, i.e. 



E'{v) 



-rjv 



1 - Qv^ 



(1-3^2)3/2' 



(21) 



where we have introduced a finite-mass dependence through the reduced mass ratio. We can see from this equation 
the LSO is found by demanding E'{v) — 0, giving viso = l/\/6. We can also see that a pole exists at Vpoie = l/VS- 
This corresponds to the last unstable circular orbit or photon ring. However, for the flux function we only have a FN 
expansion, i.e. 



■ 11 

E 

.fc=0 
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akv 



k=& 



where Fn{x) is the dominant Newtonian flux function given by 

32 



Fwix) 



(22) 



(23) 



The coefficients in the expansion of the flux function are given by 0, 0, 0, 0, 0, ^| 

1247 , 44711 8191 TT 



ao 



1, 



ai 



0, 



a2 



336 



9072 



05 



672 



6643739519 1712 7 16 tt^ 3424 ln(2) _ 16285 tt 

69854400 105~ ^ ~3 105 ' ^ 504 ' 

323105549467 2325977 1369 tt^ 39931 ln(2) 47385 ln(3) 

3178375200 4410 126 ^ 294 1568 
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2659786675197r 684877r 136967rln2 
745113600 105 105 ' 

2500861660823683 916628467 7 4242237r2 83217611 ln2 47385 ln3 

2831932303200 ^ 7858620 6804 1122660 196 

839930975040l7r 17729377r 85212837rln2 1421557rln3 

101708006400 ^ 1176 ^ 17640 784 ' 



1712 , „ , 232597 , 68487r , 916628467 , 1772937r 

Ofi = , 07 = 0, Os = , bq — , Oin = , Oil — , 25) 

105 ' ' 4410 ' 105 ' 7858620 ' 1176 ' ^ ' 

where 7 is Euler's constant. We can see that we begin to encounter logarithmic terms at fc = 6 and above. It is well 
know that terms such as these can destroy the convergence of a power series expansion. 



III. THE TEST-MASS CHEBYSHEV APPROXIMATION FLUX FUNCTION TO 5.5-PN ORDER. 

In general the Chebyshev polynomials of the first kind, T„(a;), are bound in the domain [—1, 1]. For our purposes 
we need the shifted Chebyshev polynomials, T*{x), which are bound in the domain v G [0, fi^o]. While the shifted 
polynomials naturally introduce a linear term, we found it better to forcefully introduce a linear term before expanding 
in shifted polynomials. To this extent, the first part of our analysis mirrors the steps taken in preparing the flux for 
Padeing : firstly, factor out the logarithmic terms in the PN expansion to give 



FtAv) = Fn{v) 



1 + ln 



Vis 



k=6 



11 

E 

.k=0 



C.V 



(26) 



where we have normalized the velocity in the log-terms with viso in order to minimize the effect of these terms as we 
approach the LSO and the coefficients Ck are now in general functions of viso- Secondly, we introduce a normalized 
linear term according to 



/t„(«) 



1 



^pole 



Ft„{v), 



(27) 



where Vpoie = 1/ V3, to give 



1 + ln 



V 
Viso 



11 



k=6 



E 



.k=0 



(28) 



where = c„ and f^^c^- c^^Jvpoie, fc = 1, . . . ,n. 

This is where the similarity stops with the Pade based templates. The next step is to expand both power series 
in the above equation in terms of the shifted Chebyshev polynomials. This is done by writing each monomial in the 
power series in terms of the shifted Chebyshev polynomials (see Appendix C) and substituting back into the series 
above. So, for example, the shifted Chebyshev polynomial T^{v), calculated from Equation H13() . is given by 



T*{v) = 18432^6^5 - 46O8O1;* + 6720\/6w^ - 24001;^ + 50\/6w - 1. 



(29) 



Solving this equation for , and using the pre-calculated expressions for the monomials v to v"^ in terms of the shifted 
Chebyshev polynomials, we can write 



1 

18432V6 



[Tg (w) + 10T4*(u) + 45r3*(i;) + 120r2*(z;) + 210ri*(?;) + 126To*(z;)] 



(30) 



7 



n 


5 


6 


7 


8 


9 


10 


11 


^C„ 


4.3 X 10"^ 
9.9 X 10"* 


4.7 X 10"^ 

1.8 X 10"'' 


7.8 X 10"^ 
7.6 X 10"^ 


5.4 X 10"^ 
2.6 X 10"^ 


1.7 X 10"^ 
5 X 10"^ 


1.4 X 10"^ 

8.5 X 10"^ 


2.8 X 10"^ 
7.1 X 10"* 


tot 
Tn 
tot 

c„ 


1.37 

1.3 X 10"^ 


0.94 
2.9 X 10"" 


0.47 
1.1 X 10"'' 


0.39 
3.2 X 10"^ 


0.34 
5.9 X lO"** 


0.17 

9.2 X 10"^ 


2.8 X 10"^ 
7.1 X 10"* 



TABLE L The top two lines give the truncation errors associated with each order of approximation for both the PN and 
Chebyshev flux functions. The bottom two lines give the total truncation error incurred as we reduce the order of approximation 
from 5.5 to 2 PN. 



and 



Proceeding like this for all monomials, it then allows us to write the power series in the PN expansion solely in terms 
of the shifted Chebyshev polynomials. The final expressions are 

11 11 

k=6 k=0 

401659437 742347961 345804298 345491233 

"21704658985 ° ~ 23317109996 ^ ~ 19651920683 ^ ~ 37589029195 

86276497 54224013 _ 18803603 _ 9312631 

29531863489 88082754350 221306089058 881154831507 ^ ^ 

8680152 _ 4802925 _ 785237 _ 312401 

4609576023097 ® 17752354221551 28637321914188 234989668504724 



k=0 k=0 

4343153738 2082443029 364597801 217544857 

6796318687 ° ~ 6130391137 ^ ^ 11542494373 ^ ^ 25594225319 ^ 

178665325 , 19866586 , 23168597 , 35609151 , 

T4 (u) H Tr*(v) H Tfifu) H r7*(w) (32) 

185831648387 ^ 20058344651 ' 127642232741 ^ 469435094981 ' ^ ' 

16329253 , 8570231 , 29902967 , , 1247733 

H To*(w) H T„*(i;) H TUv) H ^iV")- 

626720859540 ^ 17037955604303 ^ 3528952934787 ^ 17652943967920 ^ 

We should emphasise the fact here that although the values of the coefficients Ik are zero up to fc = 6, the Chebyshev 
expansion includes terms from fc = 0. This allows us to define the Chebyshev approximation to the gravitational wave 
flux function as 



-1 



FcAv)={^ Fn{v) 

\ Vpole J 

where we re-introduce the pole at the photon ring. 



1 + 1^(7^)E^^^"*W flKT:{v) 

^ '"^"^ k=0 n 



11 

E 

.k=0 



(33) 



IV. MODELLING THE TEST-MASS FLUX FUNCTION. 



A major reason for using Chebyshev polynomials is a process called Chebyshev economization l42j. This works as 
follows : we can always expand a power series Pn{x) of some function f{x) in terms of a Chebyshev series qn{x) such 
that 



n 



Pn{x) = E (^kX^ = E ^kTk{x) = qn{x). (34) 



fc=0 fc=0 
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r/M r/M 

FIG. 3; The error of the PN, Pade and Chebyshev approximations when compared to the numerical flux for 2 to 3.5 PN order. 



In general the series Pn{x) diverges as we move away from the point of expansion. Also, we usually find that 
I In I > |an-i|- Therefore if we wish to truncate the polynomial Pn{x) to give another polynomial rm{x), such that 
m < n, we introduce an additional truncation error of the order |a„a;"| in rjn{x)- On the other hand, the coefficients 
bk in a Chebyshev always decrease as we increase the order of approximation such that < |&„_i|. If we now try the 
same truncation with the Chebyshev series qn{x) to give another Chebyshev series Sm{x), the additional truncation 
error is guaranteed to be smaller than an equivalent truncation in the Taylor series. This is due to the fact that 
|7ri(x)| < 1 in the interval meaning that the truncation error is of the order of the truncated coefficient |6„| which 
is in general a small number. We refer the reader to Appendix A for an example of Chebyshev economization. The 
other main advantage of using Chebyshev polynomials is the fact that because they are closely related to the minimax 
polynomial, their convergence has absolutely nothing to do with the convergence of the initial Taylor series. This is 
due to the fact that all the information contained in a Taylor series is acquired from just one point in the interval, i.e. 
the point of expansion. A Chebyshev series, on the other hand, uses information from all points in the interval. It is 
these two characteristics that we believe make the Chebyshev polynomials critical for the design of any GW template. 

It has been shown in the past that for a test-mass in a circular equatorial orbit around a Schwarzschild black hole, 
that resumming the PN flux approximation using Pade approximation produces a better fit to the numerical flux. 
However, this method does suffer from certain aspects such as poles and singularities, and the fact that it may not 
perform as well as the PN flux at low orders of approximation. In this section we investigate the modelling of the 
numerical flux ^43j using Chebyshev economization. It has also been well known that the PN approximation has 
its largest error at the LSO. We can therefore estimate the truncation error involved with the PN approximation 
at various orders of magnitude. We would expect as v takes on a maximum value at the LSO, that the maximum 
truncation error should be well approximated by the first neglected term, i.e. for fT„_i(2;) the error should be smaller 
than 

eT„ = anWC^o + M'"lso)bnVuo- (35) 

On the other hand, the truncation error associated with an expansion in Ultraspherical polynomials is of the form 

ep»(a;) = c„P"(x). (36) 

For the Chebyshev approximation to the flux, we know that over the entire interval v € [0, viso] the shifted Chebyshev 
polynomials have a maximum absolute value of |T'*('(;)| < ±1. Therefore, the induced truncation error will be well 
approximated by the size of the coefficient of the truncated term, i.e. 



<5c„ < A„, 



(37) 
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r/M r/M 

FIG. 4: The error of the PN, Pade and Chebyshev approximations when compared to the numerical flux for 4 to 5.5 PN order 



due to the fact that at the LSO, the logarithmic terms are kiUed off. 

In Table I we compare the induced truncation errors in both the PN and Chebyshev expansions, assuming that 
the maximum error in both cases is at the LSO. We can see that the truncation errors involved are much smaller 
in the Chebyshev series than in the PN expansion. In fact, truncating the PN approximation from 5.5 PN to 2 PN 
order we induce a total truncation error of greater than unity. However truncation of the Chebyshev series to 2 PN 
accumulates a total error of ~ 10""^. 

In Figure (O we present the error of the PN, Pade and Chebyshev approximations when compared to a numerical 
flux at 2 to 3.5 PN order. The first thing we should point out is there in a known pathology in the PN expansion 
at 2.5 order. We can see that the PN approximation at this error performs worse than at any other order. As proof 
that the convergence of the Pade approximation is dependent on the convergence of the PN approximation, we can 
see that the Pade flux also performs worse at this order than at any other order (free from singularities). But we 
also see another downside of the Pade method in that the 2 PN approximation has a greater error at the LSO than 
the PN flux. So in this case the Pade approximation really only begins to make a difference from 3 PN onwards. 
We can see from this figure that the Chebyshev approximation tries to find an equal-ripple error curve. It achieves 
this at high values of r, but as we approach the LSO the error begins to grow. Over most of the plotted range, the 
Chebyshev approximations lie below both the PN and Pade curves. In order to do this, the Chebyshev approximation 
allows the error to grow in regions where we have good agreement in order to find a better fit elsewhere. While the 
Chebyshev approximation is not as accurate as the PN or Pade at high values of r, it has no influence on our results. 
With the test-mass systems we have chosen, the lightest system, (100, 1.4)Mq, comes into the detector bandwidth at 
about 15M, while the heaviest, (20, 1.4)M0, comes in at about 45M. We can see from the graph that in the region 
6 < M < 45, in nearly all cases, the Chebyshev approximation performs better than the PN and the Pade at all 
orders of approximation. 

In Figure Q we plot the same thing for the flux at 4 to 5.5 PN order. Here is where we see the main problem with 
the Pade approximation. While in all cases it is better than the PN approximation, there are problems. The 4.5 PN 
Pade approximation is only marginally better than the PN approximation of the same order, the 5 PN approximation 
has a singularity making this order completely unusable in the construction of any template and the 5.5 order is 
phenomenally good. Just as going to higher orders of approximation with the PN expansion does not necessarily 
guarantee us a better model, we can never be quite sure what we are going to get with the Pade approximation. 
We can see that, besides the 5.5 order where the Pade model is excellent, the Chebyshev model once again gives a 
superior fit to the numerical flux. We should point out here that the equal-ripple curve is not so obvious here as we 
are beginning to recover the original PN approximation. Even so, the Chebyshev approximation is still better than 
the original approximation at 5.5 PN order. We know that the main degradation in signal to noise ratio comes from 
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a template going out of phase at the LSO. In Figure (|3J) we plot the percentage errors in each of the approximations 
at the LSO. The PN approximation performs worst with errors of between 3.4% < ep„ < 42%, with the 2.5 PN 
approximation incurring the highest error. We can see, in fact, that the 2 PN approximation actually incurs the 
smallest error at the LSO. This is precisely the order that the Pade approximation performs worst. However, besides 
the 5.5 PN order where the error is 0.06%, the error for the Pade flux at the LSO is 2.9% < e^^^^ < 7.8%. While 
we have plotted the Pade error at 5 PN order, wc must remember that there is a singularity at this order, so the 
result is superfluous. In contrast, and a further sign that the Chebyshev approximation to the flux is related to the 
minimax polynomial, the error for the Chebyshev flux is approximately constant, 1.5% < ^cheb — 1-8%, at all orders 
of approximation. This is what we would expect as it attempts to minimize the maximum error at the LSO. 



V. RESULTS AND DISCUSSION. 



While the previous graphs show that the Chebyshev approximation provides a better fit to the numerical flux, the 
only true test is to compare templates constructed with the various approximations against some 'exact' signal. In 
order to do this we use the technique of matched filtering to see how each template performs. For our 'exact' signal 
we use a restricted PN waveform where we use the exact expression for the orbital energy function, and a numerical 
gravitational wave flux function. We define the overlap between two waveforms h{t) and s{t) as the inner product of 
the normalized waveforms denoted by 

0= , , (38) 



where the scalar product is defined by 



(/^k)=2^°°^ [h{f)r{f) + h*{f)s{f)\. (39) 

Here, an asterix denotes a complex conjugate and a tilde denotes the Fourier transform of the time domain waveform. 
Geometrically, we can think of the Fourier transform waveforms as vectors in a Hilbcrt space. As both waveforms are 
normalized, the overlap returns the cosine of the angle between them. If both waveforms match exactly, the overlap 
returns unity, but as the waveforms begin to differ, the overlap drops to zero. 

Each template is a function of a number of parameters, A'^, which in turn define the dimensionality of the search 
space. For the types of systems we are considering, as well as the extrinsic parameters Iq and 00, each template 
is defined by the individual masses of the systems, (mi,m2). Due to the short duration of the signal parameters 
such as orbital inclination, position in the sky etc. are unimportant. While we start off with a 4-d search space. 
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FIG. 6: The Fitting Factors (top) and %-error in the estimation of the chirp mass (bottom) for FN, Fade and Chebyshev 
templates as compared to an 'exact' test-mass template with parameters of (mi, 7712) — (100, 1.4)Mq. The trough in the top 
cell and the peak in the bottom cell for the Fade template at 5 FN is due to a singularity in the flux function at this order. 



= {to, (f)o,mi,m2}, we can reduce the search to the 2-d subspace of extrinsic parameters as follows : we can 
maximize over the time of arrival of the wave by simply computing the correlation of the template with the data in 
the frequency domain followed by the inverse Fourier transform. This yields the correlation of the signal with the 
data for all time-lags. The phase at time of arrival is maximized by generating two orthonormalized ternplates with 
an in-phase ^0 = and quadrature-phase (po = -k 12. The easiest way of doing this is by noticing that if ho{f) is the 
in-phase template, the quadrature template is defined by /i7r/2(/) = iho^f). The overlap maximized over the phase 
at the time of arrival is then given by 

Oma.,„ - ^J{h{^0 = 0)\s)^ + {h{^0=7r/2)\s)\ (40) 

We finally define the fitting factor J-J-' as the overlap maximized over all parameters 

^ ^ max • (^-^) 

For this study, as we are concentrating on the test-mass regime, we have chosen three systems : (100, 1.4), (50, 1.4) 
and (20,1.4) Mq. These correspond to reduced mass ratios of ?7 = 0.0136,0.0265 and 0.06 respectively. While we 
realise that the last system is pushing the test-mass regime for the ground based detectors, it still is informative as 
it shows how the fiux approximation begins to behave as we approach comparable masses. In order to calculate the 
fitting factors we used the EURO noise curve. While the waveforms are described by the two individual mass, it is 
impossible for a single ground based detector to detect these parameters. What is possible is a combination of the 
individual mass called the chirp mass. This is defined as Mc = mif'f^ . 

In Figures © to (|HJl we plot the fitting factors and percentage errors in the estimation of the chirp mass. We 
will deal with each in turn. The top cell of Figure 10 shows the fitting factors for the PN, Fade and Chebyshev 
templates for the masses (100, 1.4) Mq. In general, we will take a template to be adequate if it reaches a predetermined 
threshold. This is usually a fitting factor of about 0.97 . In this case the PN templates at 2 and 2.5 PN achieve 
fitting factors of about 0.96 and 0.87 respectively (see Appendix D for actual numbers) making them inadequate as 
templates. Only from 3 PN onwards do the templates meet the required threshold. In keeping with everything that 
we know about the PN approximation, we can see that the fitting factors follow an oscillatory pattern, i.e. the 3.5 
PN is better than both the 3 and 4 PN templates etc. We can see from the bottom cell that not only do the lower 
order approximants achieve bad fitting factors, but they also have the highest errors in the estimation of Mc- We can 
also see that the 4 PN template which performs worse than the 3.5 and 4.5 templates also has a higher parameter 
estimation error. The Fade templates achieve excellent fitting factors from 2 PN onwards, but again we can see an 
oscillatory nature in the estimation of the chirp mass. While the 2 PN template has a fitting factor of ~ 0.99 it has 
an error in the estimation of Mc almost equal to the PN template. The error in the estimation of the chirp mass 
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FIG. 7: The Fitting Factors (top) and %-error in the estimation of the chirp mass (bottom) for FN, Fade and Chebyshev 
templates as compared to an 'exact' test-mass template with parameters of (mi, 7712) — (50, 1.4)Mq . 



begins to improve as we increase the order of approximation, but we can see that the 3.5 template has a higher error 
in the chirp mass than the corresponding PN template. We should also explain that the sudden dip in the fitting 
factor and peak in the chirp mass estimation at 5 PN is due to the singularity in the flux function giving a zero fitting 
factor and infinite error in the chirp mass. At 5.5 PN order the Pade template achieves almost a perfect overlap and 
the lowest error in Ale- The Chebyshev templates, on the other hand, achieve fitting factors of almost unity at all 
orders of approximation. Not only that, but we can see that the error in estimating Mc improves, converging to a 
constant error value, which is again what we would expect as the template is near minimax. Just looking at the error 
in parameter estimation, we can see that the Chebyshev templates have achieved by 3 PN order what takes the Pade 
templates until 5.5 PN order to achieve. 

In Figure ||7J) we present the same results for the masses (50, 1.4) M©. We can see in the top cell that while the 2 PN 
template now meets the require fitting factor, the 2.5 template again falls short. Also, we notice that the oscillation 
in fitting factors is not as great as in the previous case with all PN templates achieving fitting factors of > 0.99 from 

3 PN order onwards. Once again, the Pade templates outperform the PN templates at all orders. However, we again 
see that in order to do this, it means that the Pade templates occasionally incur larger errors in the estimation of 
the chirp mass. In this case we see that the 3 PN Pade template has an error greater than the corresponding PN 
tempalate. We also see this time that the error in the estimation of Mc is only slightly better in a number of cases (2 
PN, 3.5 PN, 4 PN and 4.5 PN orders) than the equivalent PN templates. On the other hand, the Chebyshev templates 
again achieve fitting factors at all orders of approximation. We also find, possibly due to the longer waveform and the 
accumulation of phase information, that while the error at 2 PN order is roughly an order of magnitude better than 
the PN template, from 2.5 PN order onwards the best-fit template has parameters almost identical to the signal we 
are trying to fit. While hard to see from the plot, the error in the estimation of Mc is approximately 10~*% from 2.5 
to 5.5 PN order. Once again, we can see that the 2.5 PN Chebyshev template performs as well as the 5.5 PN Pade 
template and better than any of the PN templates. 

Finally, in Figure lO we present the results for the masses (20, 1.4) M©. In this case we are pushing the test-mass 
regime to the limit, therefore, we cannot be truly confident in these results. They may serve to give an idea of how 
the various approximations behave as we approach the comparable-mass case. We can see from the top cell that the 
PN templates again do not attain the necessary threshold value until the 3 PN template. What is different here to the 
other cases is that the best PN template is at the 3.5 PN order, with the fitting factors falling off afterward and only 
improving again at 5.5 PN order. In terms of estimating the chirp mass we can see that the error incurred by the PN 
templates either stays pretty much constant or matches the parameters of the signal correctly. The Pade templates 
again perform better than the PN templates in general. However, they actually achieve lower fitting factors at 3.5 and 

4 PN orders. We can see from the lower cell that they provide only marginal improvement in the estimation of the 
chirp mass. The Chebyshev templates, while not as being as good in this case as in the other two cases at 2 and 2.5 PN 
order, once again outperform both other approximations with fitting factors of > 0.99 at all orders of approximation. 
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FIG. 8: The Fitting Factors (top) and %-error in the estimation of the chirp mass (bottom) for FN, Fade and Chebyshev 
templates as compared to an 'exact' test-mass template with parameters of (mi, 7712) — (20, 1.4)Mq . 



Again, while there is an improvement in the estimation of the chirp mass at 2 PN order, the main improvement again 
comes at 2.5 PN order onwards with the templates again almost matching the parameters of the signal. While the PN 
templates provide excellent estimation of the chirp mass at 3.5, 5 and 5.5 PN orders, coupled with the corresponding 
fitting factors, the Chebyshev templates are clearly superior to the PN and Pade approximations. 

VI. CONCLUSION. 

In this study we have introduced a new template for detecting gravitational waves from compact binary systems. 
In the test-mass regime, while having an exact expression for the orbital energy function, we only have a PN approx- 
imation for the gravitational wave flux function. We have shown that one of the main problems with expansions such 
as the PN or Pade approximations is that both adhere to the Weierstrass's theorem in that we need a large number 
of terms to properly approximate a function. A better method is to use the Chebyshev Alternation theorem which 
predicts the existence of a minimax polynomial which achieves an equioscillation at at least n-\-2 points in an interval. 
While we are unable to find such a polynomial in modelling the flux function, we have shown that a member of the 
family of Ultraspherical polynomials, the Chebyshev polynomials of the flrst kind, are closely related to the minimax 
polynomial due to the fact that they achieve equioscillation at n -I- 1 points in an interval. 

We demonstrated that by using the shifted Chebyshev polynomials, we can define a new Chebyshev approximation 
to the gravitational wave flux function. A major advantage of this new flux function is that when we expand a 
series in terms of Chebyshev polynomials, each subsequent coefficient is smaller than the previous one. As the shifted 
polynomials have a maximum value of unity over the domain of interest, it means that the truncation error incurred 
by going to lower orders of approximation is proportional to the neglected coefficient, which is in general a small 
number, and is much smaller than the truncation error involved in the PN approximation. 

By graphically fitting the PN, Pade and Chebyshev approximations to the gravitational wave flux function to a 
numerical flux, we were able to show that at all orders of approximation the Chebyshev approximation provided a 
better fit than the PN approximation, and was better than the Pade approximation at all orders except the 5.5 PN 
order. The closeness of the Chebyshev flux to the minimax flux was observed due to the fact that the new flux tried 
to achieve an 'equal-ripple' error curve. One of the main results of this study is the fact the new flux function has a 
lower error at the LSO, where the templates are most likely to be out of phase with a signal, than both the PN and 
Pade approximations. 

Finally, by choosing a number of test-mass systems and a flducially 'exact' signal, we found that not only does 
the Chebyshev template always achieve higher fltting factors, but they have excellent error estimation, in most cases 
finding the parameters almost perfectly. Another of the main features of the new templates are their ability to perform 
as well at lower orders of approximation than either the PN or Pade templates at the highest order. 
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While this work is currently being extended to the comparable-mass case, it is obvious from this study that 
Chebyshev templates are superior to the PN and Pade templates in the test-mass case. We believe that they will 
prove to be an invaluable addition to any strategy that involves the detection of gravitational waves using both ground 
and space-based detectors. 
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APPENDIX A: CHEBYSHEV ECONOMIZATION. 



In this section we demonstrate the strength of using a Chebyshev economization on a power series. This process 
is sometimes also referred to as telescoping a power series. Let us say, for example, that we are looking for the best 
cubic expansion to the function cos(a;) such as 

cos(x) = a + + cx^ + O (.T^') , y x e [-1, 1] (Al) 

The first thing we would normally do is find a Taylor expansion to cos(x), i.e. 

X2 ^y>4 ry»6 rplO rrt'^ 

cos(x)«l-- + --- + --- + .. + -. (A2) 

The remainder of this series is given by (x""'"2/(n + 2)!) cos(^), and since cos{^)max = 1; any truncation gives an error 
smaller than the first neglected term. So if we truncate after x^, we introduce an error of e < |1/10!| = 2.76 x 10~^. If 
we keep going and find our cubic expansion, we find that by truncating after the x'^ term we introduce an additional 
error of e < |1/6!| = 1.4 x 10-^. 

Now, if we rewrite the Taylor series in terms of a Chebyshev series, truncated after the x^ term, we get 

cos(x) « To(.t) - [To(x) + T^ix)] + [3ro(x) + AT^ix) + T^x)] 

- 3^^[10To(x) + 15T2(x) + 6T4(x)+T6(x)] (A3) 
+ [35To(x) + 56T2(x) + 28T4(x) + 8n{x) + Tsix)] , 

which reduces to 

cos(x) « 0.7652ro(x) - 0.2298r2(x) + 4.95 x 10-^r4(x) - 4.185 x 10-^r6(x) + 1.94 x 10-^r8(x). (A4) 

Using the fact that |T'„(x)| < ±1 on the interval [—1,1], we can drop the Te{x) and Ts{x) terms, and introduce a 
maximum error of 

4.185 X 10"^ + 1.94 X 10""^ = 4.2 x 10"^ (A5) 

This is a factor of 

er 1.4 X 10-3 



ec 4.2 X 10-5 



30, (A6) 



improvement over the cubic Taylor expansion. We can see that because the Chebyshev polynomials are closely 
associated to the minimax polynomial, their convergence properties are greatly enhanced, and are clearly superior to 
a Taylor expansion of the same order. 
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APPENDIX B: T*(v) IN TERMS OF v". 

While computationally it is advantageous to be able to calculate the shifted Chebyshev polynomials using the 
recurrence relation, it is also useful to have an analytical expression of the polynomials T* (v) in terms of the monomials 
u", i.e. 

n 

T:{v)=Y.a,v''. (Bl) 

fe=0 

We present a Chebyshev table to order n = 9 below. We have written the 5 and 5.5 PN Chebyshev terms explicitly 
due to their length. 





1 


V 


2 

V 



















1 




















1 


-1 


v/24 


















2 


1 


-4\/24 


48 
















3 


-1 




-288 


192v^ 














4 


1 




960 


-1536\/6 


4608 












5 


-1 


50\/6 


-2400 


6720\/6 


-46080 


18432\/6 










6 


1 


-72^6 


5040 


-21504^/6 


248832 


-221184^6 


442368 








7 


-1 


98^/6 


-9408 


56448^/6 


-967680 


1419264\/6 


-6193152 


1769472^6 






8 


1 


128\/6 


16128 


-129024^/6 


3041280 


-6488064\/6 


46006272 


-28311552\/6 


42467328 




9 


-1 


162\/6 


-25920 


266112\/6 


-8211456 


23721984^6 


-241532928 


238878720^6 


-764411904 


169869312^6 



Ti*o(t;) = 1 - 200\/6w + 39600?;2 - 506880V6w^ + 19768320t;'* - 73801728\/6u^ + 1006387200t;^ 
- 1415577600 V^w'' + 7219445760t;^ - 3397386240\/^t;^ -|- 4076863488v^°, 

T^^{v) = -1 + 242V6V - bSOSOv^ + 906048 V6^;^ - 43490304?;'' + 202954752^6?;^ - 3542482944?;^ 

-|-6617825280\/6v^ - 47648342016?;^ + 35502686208VS?;^ - 89690996736?;^° -|- 16307453952 V^?;". 



The table is read, for example, as 

T3 (?;) = -1 + IsVev - 288v^ + 192^6?;^ 



(B2) 
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APPENDIX C: IN TERMS OF T*{v). 

In order to expand the PN series approximation as a Chebyshev series, we need to be able to write the monomials 
ti" in terms of the shifted Chebyshev polynomials according to 

n 

V- = h-'Y,'^nK{v). (CI) 
fe=0 

This then allows us to substitute the Chebyshev polynomials into the PN series. We again present the expansions in 

terms of a Chebyshev table : 





TS{v) 




TUv) 


T3*(«) 








Tr*{v) 


Ti(v) 










1 


























1 


1 






















h-zv' 


3 


4 


1 






















10 


15 


6 


1 




















35 


56 


28 


8 


1 
















b5v" 


126 


210 


120 


45 


10 


1 














bev"" 


462 


792 


495 


220 


66 


12 


1 












bjv^ 


1716 


3003 


2002 


1001 


364 


91 


14 


1 










bsv'' 


6435 


11440 


8008 


4368 


1820 


560 


120 


16 


1 










24490 


31932 


44028 


18582 


8568 


3060 


816 


153 


18 


1 








95978 


173360 


128130 


77800 


38760 


15504 


4845 


1140 


190 


20 


1 






394296 


709016 


522368 


322168 


170544 


74613 


26334 


7315 


1716 


231 


22 


1 



where the coefBcients 6„ are given by 

bo = 1, bi = V24, b2 = 48, 63 = 192\/6, 

64 = 4608, 65 = 18432^6, be = 442368, 67 = 1769472^/6, (C2) 

68 = 42467328, 69 = 169869312VS, 610 = 4076863488, 611 = 16307453952V^, 

and again, for example, the table is read as 

bav"^ = WT*{v) + 15T*{v) + 6T;{v) + T;{v). (C3) 
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APPENDIX D: FITTING FACTORS AND PARAMETER ESTIMATION 

We present the values of the fitting factors and the error in the estimation of the Chirp Mass in the foUowing tables. 
The labels T, P and C stand for the PN, Padc and Chebyshev results. Row one of each order of approximation denotes 
the fitting factor, in row two are the individual masses associated with each fitting factor and row three corresponds 
to the percentage error in the estimation of the Chirp Mass. The blank entries at 5 PN order for the Pade templates 
is due to a singularity at this order of approximation in the gravitational wave flux function. 





(100, 1A)Mq 


(50, 1.4)M0 


(20, 1.4)M0 




T 


P C 


T P 


C 


T P 


C 


/I 

4 


0.9614 


0.9895 0.9998 


0.974 0.9863 


0.9949 


0.9436 0.9815 


0.9894 




(98.83, 1.4) 


(102.13, 1.39) (99.54, 1.41) 


(49.03,1.41) (51.72,1.37) (50.08,1.4) 


(19.16,1.44) (20.86,1.36) (20.05,1.39) 




0.473 


0.422 0.24 


0.37 0.083 


0.065 


0.122 0.035 


0.023 


r 



0.866 


0.995 0.9995 


0.602 0.996 


0.998 


0.649 0.9864 


0.9929 




(116.19,1.3) 


(100.56, 1.4) (100.03, 1.4) 


(56.15, 1.35) (50.56, 1.39) 


(50, 1.4) 


(23.65, 1.25) (20.21, 1.39) 


(20, 1.4) 




1.62 


0.225 0.012 


2.57 0.025 





0.072 0.01 





D 


0.973 


0.9969 0.9992 


0.9899 0.9985 


0.9986 


0.9823 0.9868 


0.9985 




(97.52,1.41) 


(100.96, 1.39) (100.01, 1.4) 


(47.98,1.44) (50.48,1.39) (49.99,1.4) 


(19.05, 1.45) (20.19, 1.39) 


(20, 1.4) 




0.585 


0.043 0.004 


0.002 0.039 





0.045 0.03 





'7 
/ 


0.9975 


0.998 0.9993 


0.9955 0.9989 


0.9984 


0.9958 0.9917 


0.9995 




(100.79,1.39) (100.55,1.39) (100.01,1.4) 


(50.41, 1.39) (50.21, 1.39) 


(50, 1.4) 


(20, 1.4) (20.22, 1.39) 


(20, 1.4) 




0.111 


0.207 0.004 


0.095 0.211 





0.011 





8 


0.9843 


0.9969 0.9993 


0.9914 0.9979 


0.9984 


0.9868 0.9843 


0.9998 




(103.91,1.36) (100.99,1.39) (100.01,1.4) 


(51.95, 1.36) (50.49, 1.39) 


(50, 1.4) 


(20.6, 1.37) (20.2, 1.39) 


(20, 1.4) 




0.187 


0.031 0.004 


0.173 0.031 





0.052 0.009 





9 


0.999 


0.9967 0.9993 


0.9984 0.9981 


0.9985 


0.9821 0.9891 


0.9998 




(98.25, 1.42) 


(100.92,1.39) (100.01,1.4) 


(49.1,1.42) (50.47,1.39) 


(50, 1.4) 


(19.81, 1.41) (20.19, 1.39) 


(20, 1.4) 




0.136 


0.059 0.004 


0.107 0.047 





0.023 0.03 





10 


0.996 


0.9993 


0.9957 


0.9985 


0.982 


0.9998 




(100.88, 1.39) 


(-,-) (100.01,1.4) 


(50.45, 1.39) (-, -) 


(50, 1.4) 


(20,1.4) (-,-) 


(20, 1.4) 




0.075 


0.004 


0.063 











11 


0.9969 


0.9999 0.9993 


0.9927 0.9995 


0.9919 


0.9999 1.0 


0.9998 




(100.78, 1.39) 


(100, 1.4) (100.01, 1.4) 


(50.43, 1.39) (50, 1.4) 


(50, 1.4) 


(20, 1.4) (20, 1.4) 


(20, 1.4) 




0.115 


0.004 0.004 


0.079 












